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The  purpose  of  this  research  was  to  investigate  the  use  of  finite  clement  tech- ) 
niques  for  solving  convective  heat  transfer  problems  which  involve  phase  chance J 
The  specific  application  was  to  permafrost  mdoetino.  This  has  been  cccm*!  i$he<| 
at  several  levels,  ranging  from  very  theoretical  error  estimates  and  develop-  l 
ment  of  special  numerical  techniques  to  handle  phase  change,  to  application  of  ' 
finite  element  techniques  to  specific  permafrost  problems.  <-  1 
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The  purpose  of  this  project  was  to  Invest lgsts  the  use  of  finite  eleaent 
techniques  for  solving  convective  heat  transfer  p rob leas  which  Involve  phase 
change.  The  specific  application  was  to  peraafrost  aodellng.  This  project  has 
been  accoapllshed  at  several  levels,  ranging  froa  very  theoretical  error 
estiaates  and  developaent  of  special  nuaerlcal  techniques  to  handle  phase 
change,  to  application  of  finite  eleaent  techniques  to  specific  peraafrost 
problems.  Our  research  haa  focused  on  the  following  areas: 

(1)  Analysis  of  convective  effects  In  convection-diffusion  processes; 

(2)  Use  and  analysis  of  upwinded  finite  elements  to  stabilise 
convection-diffusion  models  in  the  presence  of  layers  and  fronts; 

(1)  Error  analysis  of  aovlng  noshes  to  track  a  tine-dependent  freese  front; 

(4)  Special  techniques  to  track  the  freese  front  and  to  accurately  calculate 
the  energy  contribution  due  to  the  latent  heat  of  fusion  during  phase 
change. 

(5)  Error  estiaates  for  the  finite  eleaent  alternating-direction  procedures 
developed  here; 

(()  Basic  questions  relating  to  the  maerical  solution  of  general  heat  transfer 
processes. 

D1SCPSS10K 

An  eigenvalue  analysis  (1,4)  was  used  to  study  the  theoretical  properties  of 
general  convection-diffusion  problems.  The  specific  structure  of  the 
sealdlscrete  systems  obtained  by  alternative  finite  eleaent  methods  are  analysed 
using  Cershgorin  theory  and  the  theory  of  oscillation  matrices.  This  provided 
an  interesting  approach  for  analysing  the  effects  of  higher-order  elements, 
upwlnding  and  aesh  refinement.  This  provided  a  better  understanding  of 
oscillatory  and  dissipative  effects  artificially  introduced  by  particular 
numerical  models.  We  have  also  Investigated  freesing  in  the  presence  of 
internal  convection  (16,22,23).  We  have  found  that  one  can  control 
oscillations,  in  this  case,  to  some  degree,  by  varying  the  mesh  else  and  element 
placement.  In  some  related  studies  we  have  examined  the  effect  of  different 
boundary  conditions  on  the  developaent  of  oscillations  due  to  a  a  discontinuity 
in  the  initial  data  and  reflection  of  a  front  at  a  boundary  | 9  ).  Roth 
dissipative  backward-type  Integration  schemes  and  nondisslpatlve  central  schemes 
have  been  used  in  a  comparative  nuaerlcal  studies.  The  algorithms  are  based  on 
variable  order  -  variable  step  Integration  techniques. 
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One  of  Ch«  major  difficult!**  In  dealing  with  transport  problems  in  which 
convection  plays  an  important  role  is  the  effective  treatment  of  sharp  fronts  in 
the  solution  profiles.  Some  techniques,  although  accurate,  lead  to  oscillatory 
errors  while  others  are  over-dissipative.  In  our  investigations  of  these 
problems  we  have  developed  a  class  of  upwind  finite  element  methods  in  which  the 
amount  of  upwindlng  is  controlled  adaptively  during  the  solution  procedure  (  1  ). 

The  eigenvalue  analysis  In  (  4  )  forms  the  basis  of  this  new  method  for 
variable  upwindlng  in  finite  element  solution  of  convection-dominated  flows 
(24 I.  The  approach  seems  quite  general  and  can  be  used  for  non-uniform  meshes, 
a  point  that  has  not  been  generally  true  of  other  upwind  techniques.  Ue  have 
extended  the  technique  to  nonlinear  convection-diffusion  problems  using  local 
linearisation  In  a  given  tlmestep  to  select  the  local  dissipation  in  the 
upwindlng  strategy  (  9  |.  Earlier,  good  results  were  obtained  with  linear  model 
problems.  Subsequently,  we  have  applied  the  technique  successfully  for  a 
multicomponent-multiphase  coupled  flow,  with  a  sharp  front  *211.  The  scheme  is 
also  being  used  in  conjunction  with  adaptive  mesh  refinement  for  problems  with 
sharp  layers  such  as  those  encountered  In  the  application  of  Interest  in  our 
present  research  (17).  Some  related  results  were  presented  as  an  invited  paper 
at  the  SES  meeting  in  December,  19A0. 

doit*  a  bit  of  progress  has  been  made  in  the  development  and  analysis  of 
moving  mesh  techniques  for  solving  heat  transfer  problems  both  with  strong 
convection  and  with  mild  convective  effects.  The  problem  may  be  viewed  as  an 
optimal  control  problem  In  which  both  the  solution  values  and  the  nodal 
coordinates  for  the  nesh  are  to  be  determined  at  each  tlmestep. 

A  fundamental  question  in  the  analysis  of  the  permafrost  problem  and 
similar  problems  involving  a  change  of  phase  is  the  treatment  of  the  free 
surface.  A  significant  part  of  the  research  effort  has  been  directed  towards 
the  tnplementatlon  of  a  practical  front  tracking  technique  and  of  accurate 
latent  heat  calculations  for  two-dimensional  thawing  around  a  pipe.  Numerical 
calculations  are  very  sensitive  to  the  latent  heat  contributions,  so  one  must 
balance  the  desired  accuracy  with  the  computational  expense.  The  front  location 
and  latent  heat  contribution  may  oscillate  when  standard  tine  Integration 
techniques  are  employed.  To  circumvent  this  problem  a  modified 
predictor-corrector  strategy  has  been  developed  and  implemented  in  the  two- 
dimensional  permafrost  code  (20).  A  series  of  numerical  experiments  have  been 
conducted  to  evaluate  and  compare  these  methods  with  more  standard  techniques. 
The  difficulties  involved  with  latent  heat  are  more  apparent  in  problems 
Involving  several  materials  and  irregular  geometries.  Me  Investigated  the 
effects  of  using  various  sizes  of  finite  elements,  different  orders  of  finite 
elements,  and  different  numbers  of  elements  within  a  given  material  property 
upon  the  numerical  calculations  (23,26].  This  leads  to  some  general  guidelines 
for  setting  up  finite  element  grids  in  irregular  geometries  when  latent  heat 
must  he  dealt  with  (25,26). 

Error  estimates  for  nonlinear  parabolic  p.d.e.'s  using  an  alternating- 
direction  technique  have  been  obtained  for  several  numerical  algorithms.  These 
algorithms  have  been  compared  In  terms  of  execution  time,  solution  accuracy  and 
storage  requirements  for  simple  two-dimensional  model  problems  (7,14,5].  This 
work  generalizes  the  alternating-direction  method  to  general  curved  geometries 
using  isoparametric  finite  elements.  This  technique  allows  two-  and  three- 
dimensional  problems  to  be  solved  as  a  series  of  one-dimensional  problems 
(1 3,14,6].  These  techniques  require  considerably  less  computer  time  and  storage 


than  standard  finite  difference  or  finite  element  methods.  Ue  have  Investigated 
the  use  of  reduced  Integration  to  further  Improve  computational  efficiency. 
This  lead  to  a  criteria  for  the  lowest  order  stable  nunericr'  Integration  scheme 
|12,15I.  This  criteria  depends  on  the  differential  equation  as  well  as  the 
particular  boundary  conditions. 

Finally,  part  of  our  research  efforts  were  directed  at  a  variety  of  basic 
questions  dealing  with  numerical  solution  of  complex  heat  transfer  processes, 
such  as  in  what  form  should  the  t lae-dependent  finite  element  equations  be 
written  to  Insure  that  the  proper  Interface  condition  Is  Implemented  at  material 
Interfaces  as  well  as  at  Ice-water  interfaces  (19),  and  what  options  may  one  use 
to  Implement  a  convection  boundary  condition  and  still  maintain  numerical 
stability.  Ocher  related  work  has  dealt  with  the  selection  of  appropriate  flux 
conditions  and  treatment  of  the  capacity  coefficient  In  the  governing  equation. 
One  approach  to  modeling  heat  transfer  In  permafrost  Is  to  replace  the  Infinite 
boundary  by  a  finite  boundary  and  to  apply  the  boundary  condition  at  Infinity  to 
this  finite  boundary.  This  can  lead  to  serious  errors  if  the  finite  boundary  Is 
too  close  to  the  source.  Ue  have  Investigated  the  use  of  a  mixed-type  boundary 
condition  at  the  finite  remote  boundaries  and  have  found  that  it  produces  a 
great  increase  In  accuracy  over  applying  either  a  fixed  temperature  or  a  no  heat 
flow  condition.  This  will  allow  more  accurate  modeling  over  much  smaller 
physical  domains. 

Fart  of  our  work  has  dealt  with  fictitious  oscillations  that  arise  In 
finite  element  modeling  of  transport  problems.  These  oscillations  occur  in  the 
convection  terns.  Both  aspects  have  been  studied  and  results  presented  at  a 
recent  meeting  on  numerical  methods  for  heat  transfer  (17).  With  respect  to  the 
front  oscillations,  these  difficulties  have  been  overcome  by  means  of  a  new 
modified  predictor-corrector  strategy  (20).  Problems  with  convection  are 
treated  using  a  variable  upwind  finite  element  technique  which  is  utilised  and 
demonstrated  to  be  effective  in  suppressing  oscillations,  while  not  inducing 
excessive  dissipation  and  degrading  of  the  solution  f  9  I •  Another  source  of 
numerical  oscillations  Is  when  boundary  temperatures  are  changing  faster  than 
heat  can  be  conducted  through  the  soil.  Ue  have  investigated  ways  In  which  one 
can  control  these  oscillations  by  varying  clement  site,  by  placing  small 
elements  near  the  boundary  and  by  varying  the  range  of  temperatures  over  which 
latent  heat  Is  liberated  (23).  When  transient  problems  involve  sharp  solution 
gradients  that  propagate  through  a  material  interface  or  through  an  abrupt 
change  In  mesh  site,  similar  difficulties  arise.  Ue  have  considered  some 
prototype  transient  problems  and  the  effect  of  different  boundary  conditions  on 
developswnt  of  oscillations  due  to  phase  errors  on  contact  and  reflection  of  the 
solution  (  9  ).  The  effect  of  different  ODF.  integration  techniques  on 
suppression  of  oscillatory  errors  was  also  studied  [3,10l. 

In  sissmary,  many  Important  theoretical  and  practical  questions  related  to 
efficient  and  accurate  permafrost  modeling  have  been  addressed  during  the  course 
of  this  three-year  contract.  The  results  of  this  project  have  been  widely 
disseminated  in  technical  Journal*  and  reports  and  at  technical  meetings.  The 
Investigators  have  been  very  productive  and  have  answered  many  questions  related 
to  general  convective  heat  transfer  with  phase  change. 
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